rm(list = ls())

## Packages

library(tidyverse)
library(stargazer)
library(lfe)
library(broom)

source('code/function_tidy_felm.R')

## Def helper function

str_filter <- function (string, pattern, neg = F) 
{
  if (neg == F) {
    string %>% .[str_detect(., pattern)]
  }
  else {
    string %>% .[!str_detect(., pattern)]
  }
}

## Get data

df <- read_rds("data/Data_Clean.rds") %>% 
  mutate(year = as.numeric(year)) %>% 
  mutate(unit = coalesce(ags, vb_key)) %>% 
  mutate(treat_weakly = ifelse(treat_categ == 'schwach', 1, 0),
         treat_heavy = ifelse(treat_categ %in% c('schwer', 
                                                 'sehr schwer'), 1, 0)) %>% 
  mutate(treat_weakly_post = treat_weakly * post,
         treat_heavy_post = treat_heavy * post) %>% 
  mutate(year_land = paste0(year, '_', land))

## Define treatment variables 

outcomevars <- c("turnout", "cdu_csu", "spd", "fdp", "afd", "greens", 
                 "left")
treatvars <- c("treat_weakly", 
               "treat_heavy",
               'treat_wounded_or_dead')
covars <- df %>% 
  colnames() %>% str_filter("county")

#### Function to generate differenced DFs ####

do_diff <- function(df, y, diff_years = c('2017', '2021')) {
  
  ## Get relevant variables
  
  df_out <- df %>% 
    arrange(unit, year) %>% 
    dplyr::select(unit, year, matches('treat'), !!y, county_id,
                  land, one_of(covars)) %>% 
    filter(year %in% diff_years)
  
  ## Declare outcome
  
  df_out[, 'y'] <- df_out[, y]
  
  ## Difference

  df_out <- df_out %>% 
    group_by(unit) %>% 
    summarise(y_diff = y[2] - y[1],
              county_age = county_age[2] - county_age[1], 
              county_unem = county_unem[2] - county_unem[1], 
              county_pop_tot = county_pop_tot[2] - county_pop_tot[1], 
              county_pop_foreign = county_pop_foreign[2] - county_pop_foreign[1], 
              county_foreign_share = county_foreign_share[2] - county_foreign_share[1], 
              county_abitur_share = county_abitur_share[2] - county_abitur_share[1],
              treat_weakly = first(treat_weakly),
              treat_heavy = first(treat_heavy),
              county_id = first(county_id),
              land = first(land),
              treat_wounded_or_dead = first(treat_wounded_or_dead)) %>%
    ungroup() %>% 
    dplyr::select(unit, y_diff, matches('treat'), county_id, land,
                  matches('county')) %>% 
    dplyr::rename(!!y := y_diff) %>% 
    filter(!is.na(!!y)) %>% 
    distinct()
  
  ## Return 
  
  df_out
}

#### Table 1 - main results ####

diff_df <- do_diff(df, 'greens') %>% 
  mutate(greens = greens * 100) %>% 
  mutate(treat_either = ifelse(treat_weakly == 1 | treat_heavy == 1,
                               1, 0)) 

## Models

m1 <- felm(greens ~ treat_weakly + treat_heavy | 0| 0| county_id,
           data = diff_df)
m2 <- felm(greens ~ treat_either | 0| 0| county_id,
           data = diff_df)
m3 <- felm(greens ~ treat_wounded_or_dead | 0| 0| county_id,
           data = diff_df)
m4 <- felm(greens ~ treat_weakly + treat_heavy | land| 0| county_id,
           data = diff_df)
m5 <- felm(greens ~ treat_either | land| 0| county_id,
           data = diff_df)
m6 <- felm(greens ~ treat_wounded_or_dead | land| 0| county_id,
           data = diff_df)

## List of models 

mlist <- list(m1,m4,m2,m5,m3,m6)

## Mean of DV

dv_means <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_mean) %>% 
  round(2)

## SD of DV

dv_sds <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_sd) %>% 
  round(2)

## Table 

stargazer::stargazer(list(m1,m4,m2,m5,m3,m6), keep = 'treat',
                     covariate.labels = c('Flooding: weakly affected',
                                          'Flooding: highly affected',
                                          "Flooding: any",
                                          'Any wounded or dead'),
                     style = 'ajps',
                     dep.var.labels = 'DV: Change in Green party vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'),
                     add.lines = list(c('DV mean', dv_means),
                                      c('DV SD', dv_sds)))

#### Table A.5 - additional covariates ####

diff_df <- do_diff(df, 'greens') %>% 
  mutate(greens = greens * 100) %>% 
  mutate(treat_either = ifelse(treat_weakly == 1 | treat_heavy == 1,
                               1, 0)) 

## Models

m1 <- felm(greens ~ treat_weakly + treat_heavy +
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share| 0| 0| county_id,
           data = diff_df)
m2 <- felm(greens ~ treat_either+
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share | 0| 0| county_id,
           data = diff_df)
m3 <- felm(greens ~ treat_wounded_or_dead+
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share | 0| 0| county_id,
           data = diff_df)
m4 <- felm(greens ~ treat_weakly + treat_heavy+
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share | land| 0| county_id,
           data = diff_df)
m5 <- felm(greens ~ treat_either+
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share | land| 0| county_id,
           data = diff_df)
m6 <- felm(greens ~ treat_wounded_or_dead+
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share | land| 0| county_id,
           data = diff_df)

## ## ##

mlist <- list(m1,m4,m2,m5,m3,m6)

## Mean of DV

dv_means <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_mean) %>% 
  round(2)

## SD of DV

dv_sds <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_sd) %>% 
  round(2)

## Table 

stargazer::stargazer(list(m1,m4,m2,m5,m3,m6), keep = 'treat',
                     covariate.labels = c('Flooding: weakly affected',
                                          'Flooding: highly affected',
                                          "Flooding: any",
                                          'Any wounded or dead'),
                     style = 'ajps',
                     dep.var.labels = 'DV: Change in Green party vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'),
                     add.lines = list(c('DV mean', dv_means),
                                      c('DV SD', dv_sds)))

#### Table A.8 - results by state ####

diff_df <- do_diff(df, 'greens') %>% 
  mutate(greens = greens * 100) %>% 
  mutate(treat_either = ifelse(treat_weakly == 1 | treat_heavy == 1,
                               1, 0)) 

diff_df$land %>% table()

m1 <- felm(greens ~ treat_weakly + treat_heavy | 0| 0| county_id,
           data = diff_df)
m2 <- felm(greens ~ treat_weakly + treat_heavy | 0| 0| county_id,
           data = diff_df%>% 
             filter(land == 'NRW'))
m3 <- felm(greens ~ treat_weakly + treat_heavy | 0| 0| county_id,
           data = diff_df %>% 
             filter(land == 'RP'))
m4 <- felm(greens ~ treat_wounded_or_dead | 0| 0| county_id,
           data = diff_df)
m5 <- felm(greens ~ treat_wounded_or_dead | 0| 0| county_id,
           data = diff_df%>% 
             filter(land == 'NRW'))
m6 <- felm(greens ~ treat_wounded_or_dead | 0| 0| county_id,
           data = diff_df %>% 
             filter(land == 'RP'))

## List of models

mlist <- list(m1,m4,m2,m5,m3,m6)

## Mean of DV

dv_means <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_mean) %>% 
  round(2)

## SD of DV

dv_sds <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_sd) %>% 
  round(2)

## Table 

stargazer::stargazer(list(m1,m4,m2,m5,m3,m6), keep = 'treat',
                     covariate.labels = c('Weakly affected',
                                          'Highly affected',
                                          'Any wounded or dead'),
                     style = 'ajps',
                     dep.var.labels = 'DV: Change in Green party vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'),
                     add.lines = list(c('DV mean', dv_means),
                                      c('DV SD', dv_sds)))

#### Table A.9 - AfD  ####

diff_df_afd <- do_diff(df = df, y = 'afd', 
                       diff_years = c('2017', '2021')) %>% 
  mutate(afd = afd * 100) %>% 
  mutate(treat_either = ifelse(treat_weakly == 1 | treat_heavy == 1,
                               1, 0)) 

## Models 

m1 <- felm(afd ~ treat_weakly + treat_heavy | 0| 0| county_id,
           data = diff_df_afd)
m2 <- felm(afd ~ treat_either | 0| 0| county_id,
           data = diff_df_afd)
m3 <- felm(afd ~ treat_wounded_or_dead | 0| 0| county_id,
           data = diff_df_afd)
m4 <- felm(afd ~ treat_weakly + treat_heavy | land| 0| county_id,
           data = diff_df_afd)
m5 <- felm(afd ~ treat_either | land| 0| county_id,
           data = diff_df_afd)
m6 <- felm(afd ~ treat_wounded_or_dead | land| 0| county_id,
           data = diff_df_afd)

## List of models

mlist <- list(m1,m4,m2,m5,m3,m6)

## Mean of DV

dv_means <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_mean) %>% 
  round(2)

## SD of DV

dv_sds <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_sd) %>% 
  round(2)

## Table 

stargazer::stargazer(list(m1,m4,m2,m5,m3,m6), keep = 'treat',
                     covariate.labels = c('Flooding: weakly affected',
                                          'Flooding: highly affected',
                                          "Flooding: any",
                                          'Any wounded or dead'),
                     style = 'ajps',
                     dep.var.labels = 'DV: Change in Green party vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'),
                     add.lines = list(c('DV mean', dv_means),
                                      c('DV SD', dv_sds)))

#### Table A.10 - AfD w/ additional covars ####

m1 <- felm(afd ~ treat_weakly + treat_heavy+
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share | 0| 0| county_id,
           data = diff_df_afd)
m2 <- felm(afd ~ treat_either+
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share | 0| 0| county_id,
           data = diff_df_afd)
m3 <- felm(afd ~ treat_wounded_or_dead+
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share | 0| 0| county_id,
           data = diff_df_afd)
m4 <- felm(afd ~ treat_weakly + treat_heavy+
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share | land| 0| county_id,
           data = diff_df_afd)
m5 <- felm(afd ~ treat_either+
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share | land| 0| county_id,
           data = diff_df_afd)
m6 <- felm(afd ~ treat_wounded_or_dead+
             county_age+county_unem+county_pop_tot+
             county_foreign_share+
             county_abitur_share | land| 0| county_id,
           data = diff_df_afd)

## ## ##

mlist <- list(m1,m4,m2,m5,m3,m6)

## Mean of DV

dv_means <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_mean) %>% 
  round(2)

## SD of DV

dv_sds <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_sd) %>% 
  round(2)

## Table 

stargazer::stargazer(list(m1,m4,m2,m5,m3,m6), keep = 'treat',
                     covariate.labels = c('Flooding: weakly affected',
                                          'Flooding: highly affected',
                                          "Flooding: any",
                                          'Any wounded or dead'),
                     style = 'ajps',
                     dep.var.labels = 'DV: Change in Green party vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'),
                     add.lines = list(c('DV mean', dv_means),
                                      c('DV SD', dv_sds)))

#### Table A.11 - all periods ####

df <- df %>% 
  mutate(treat_either = ifelse(treat_weakly == 1 | treat_heavy == 1,
                               1, 0)) %>% 
  mutate(treat_either_post = treat_either*post,
         treat_wounded_or_dead_post = treat_wounded_or_dead*post)

## Models 

m1 <- felm(I(greens*100) ~ treat_weakly_post + treat_heavy_post | 
             county_id +year | 0 | county_id, data = df)
m2 <- felm(I(greens*100) ~ treat_weakly_post + treat_heavy_post | 
             county_id +year_land | 0 | county_id, data = df)
m3 <- felm(I(greens*100) ~ treat_either_post | 
             county_id +year | 0 | county_id, data = df)
m4 <- felm(I(greens*100) ~ treat_either_post | 
             county_id +year_land | 0 | county_id, data = df)
m5 <- felm(I(greens*100) ~ treat_wounded_or_dead_post | 
             county_id +year | 0 | county_id, data = df)
m6 <- felm(I(greens*100) ~ treat_wounded_or_dead_post | 
             county_id +year_land | 0 | county_id, data = df)

mlist <- list(m1,m2,m3,m4,m5,m6)

## Mean of DV

dv_means <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_mean) %>% 
  round(3)

## SD of DV

dv_sds <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_sd) %>% 
  round(3)

stargazer::stargazer(mlist, keep = 'post',
                     covariate.labels = c('Flooding: weakly affected',
                                          'Flooding: highly affected',
                                          "Flooding: any",
                                          'Any wounded or dead'),
                     style = 'ajps',
                     dep.var.labels = 'DV: Green party vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'), type = 'latex')

#### Table A.12 - placebo ####

diff_df <- do_diff(df = df, y = 'greens', 
                   diff_years = c('2013', '2017')) %>% 
  mutate(greens = greens * 100) %>% 
  mutate(treat_either = ifelse(treat_weakly == 1 | treat_heavy == 1,
                               1, 0)) 

##

m1 <- felm(greens ~ treat_weakly + treat_heavy | 0| 0| county_id,
           data = diff_df)
m2 <- felm(greens ~ treat_either | 0| 0| county_id,
           data = diff_df)
m3 <- felm(greens ~ treat_wounded_or_dead | 0| 0| county_id,
           data = diff_df)
m4 <- felm(greens ~ treat_weakly + treat_heavy | land| 0| county_id,
           data = diff_df)
m5 <- felm(greens ~ treat_either | land| 0| county_id,
           data = diff_df)
m6 <- felm(greens ~ treat_wounded_or_dead | land| 0| county_id,
           data = diff_df)

mlist <- list(m1,m4,m2,m5,m3,m6)

## Mean of DV

dv_means <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_mean) %>% 
  round(3)

## SD of DV

dv_sds <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_sd) %>% 
  round(3)

## Table 

stargazer::stargazer(list(m1,m4,m2,m5,m3,m6), keep = 'treat',
                     covariate.labels = c('Flooding: weakly affected',
                                          'Flooding: highly affected',
                                          "Flooding: any",
                                          'Any wounded or dead'),
                     style = 'ajps',
                     dep.var.labels = 'DV: Change in Green party vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'),
                     add.lines = list(c('State FE', rep(c('No', 'Ye'), 3)),
                                      c('DV mean', dv_means),
                                      c('DV SD', dv_sds)))

#### Table A.14 - municipal-level treatment ####

## Note that this table requires data that we cannot make available publicly

diff_df <- do_diff(df = df, y = 'greens', 
                   diff_years = c('2017', '2021')) %>% 
  mutate(greens = greens * 100)

## Load treatment data

st <- read_rds('data_private/treatment_muni_level.rds')

## Merge

diff_df <- diff_df %>% 
  left_join(st, by = c('unit' ='muni_id'))

## Models

m1 <- felm(greens ~ gem_any_flooded | 0| 0| county_id,
           data = diff_df)
m2 <- felm(greens ~ gem_any_flooded | land| 0| county_id,
           data = diff_df)
m3 <- felm(greens ~ gem_any_flooded | county_id | 0| county_id,
           data = diff_df)

##

mlist <- list(m1, m2, m3)

## MEan of DV

dv_means <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_mean) %>% 
  round(2)

## SD of DV

dv_sds <- mlist %>% 
  lapply(tidy_felm) %>% 
  reduce(rbind) %>% 
  filter(term == '(Intercept)') %>% 
  pull(dv_sd) %>% 
  round(2)

## Table

stargazer::stargazer(list(m1,m2,m3), keep = 'gem',
                     covariate.labels = 'Flooding in municipality',
                     style = 'ajps',
                     dep.var.labels = 'DV: Change in Green party vote share (p.p.)',
                     model.numbers = F, 
                     keep.stat = c('n', 'rsq'),
                     add.lines = list(c('DV mean', dv_means),
                                      c('DV SD', dv_sds)),
                     type = 'latex')


